# Author: Ivan Canay | Northwestern University | Nov 2019
# Copyright (c) 2019 Northwestern University. All rights reserved.
#-------------------------------------------------------------------
# X-plain: 	This file has several function in R that are used in the
#          	joint project with Azeem M. Shaikh and Andres Santos.
#			"The wild bootstrap with a small number of large clusters"
#-------------------------------------------------------------------
require("Matrix");
require("mvtnorm")
require("ggm")
require("dummies")
require(sandwich, quietly = TRUE)
require(lmtest, quietly = TRUE)
options(warn=-1)

#-------------------------------------------------------------------
gen.data.one <-function(q,n,design=1,beta=1){
#-------------------------------------------------------------------
# X-plain: Generates data for the simulations
#-------------------------------------------------------------------
# INPUTS: - q: the number of clusters
#         - n: the reference sample size within cluster (n=0 for hete)
#         - design: sets the dgp
#------------------------------------------------------------------
# RETURNS: df = (cluster, unit, Yij, Zij)
#   Clust Unit Yij       Zij      
#1    1    1  4.0651708 0.7475937 ...
#2    1    2 -0.7227345 0.2987556 ...
#3    1    3  2.4296410 0.5015247 ...
#4    1    4  0.3724746 0.4901809 ...
#5    2    1 -0.6716404 0.2786336 ...
#6    2    2  1.6568098 0.1311757 ...
#------------------------------------------------------------------
# Setup sample sizes per cluster
half.1= floor(q/2);
half.2= q-floor(q/2);
if (n>0){
    N=n*q;
    n_j = c(rep(n,q));
} else {
    n_j = c(rep(50,each=half.1),rep(300,each=half.2));    
    N = sum(n_j);
}
# Cluster indicator (n>0 for n_j=n_j' - n=0 for n_j != n_j')
if (n>0){
    cluster = rep(c(1:q),each=n);
    unit    = rep(c(1:n),q);
} else {
    cluster=rep(c(1:q),n_j)
    unit.1    = rep(c(1:50),half.1);
    unit.2    = rep(c(1:300),half.2);
    unit    = c(unit.1,unit.2);
}
# Generate regressors
    X_j  = rep(rnorm(q),n_j);
    X_ij = rnorm(N);

#Generate Error terms
    delta_j     = rep(rnorm(q),n_j);
    epsilon_ij  = rnorm(N);

# Final random variables by design
if (design==1){
    Z    = X_j + X_ij;
    U    = (Z^2)*(delta_j + epsilon_ij);
} else if (design==2){
    Z    = sqrt(cluster)*(X_j + X_ij);
    U    = (Z^2)*(delta_j + epsilon_ij);
}
## Designs below are the same as 1 and 2 but with different U
if (design==11){
    Z    = X_j + X_ij;
    U    = delta_j + (epsilon_ij*X_ij^2 + X_ij^2-1);
} else if (design==12){
    Z    = sqrt(cluster)*(X_j + X_ij);
    U    = delta_j + (epsilon_ij*X_ij^2 + X_ij^2-1);
}
# Generate Outcome
    Y    = 1 + beta*Z + U;

# Return the data frame
    df  = data.frame(cluster,unit,Y,Z);
    return(df);
}

#-------------------------------------------------------------------
gen.data.k <-function(q,n,corr,design=3,beta1=1){
#-------------------------------------------------------------------
# X-plain: Generates data for the simulations
#-------------------------------------------------------------------
# INPUTS: - design: ranges from 3 to 4.
#         - Z: covariates generated by gen.Z 
#         - q: the number of clusters
#         - n: the reference sample size within cluster 
#------------------------------------------------------------------
# RETURNS: df = (cluster, unit, Yij, Zij)
#   Clust Unit Yij       Zij      
#1    1    1  4.0651708 0.7475937 ...
#2    1    2 -0.7227345 0.2987556 ...
#3    1    3  2.4296410 0.5015247 ...
#4    1    4  0.3724746 0.4901809 ...
#5    2    1 -0.6716404 0.2786336 ...
#6    2    2  1.6568098 0.1311757 ...
#------------------------------------------------------------------
# Setup sample sizes per cluster
half.1= floor(q/2);
half.2= q-floor(q/2);
if (n>0){
    N=n*q;
    n_j = c(rep(n,q));
} else {
    n_j = c(rep(50,each=half.1),rep(300,each=half.2));    
    N = sum(n_j);
}
# Cluster indicator (n>0 for n_j=n_j' - n=0 for n_j != n_j')
if (n>0){
    cluster = rep(c(1:q),each=n);
    unit    = rep(c(1:n),q);
} else {
    cluster=rep(c(1:q),n_j)
    unit.1    = rep(c(1:50),half.1);
    unit.2    = rep(c(1:300),half.2);
    unit    = c(unit.1,unit.2);
}
#Generate Error terms
    delta_j     = rep(rnorm(q),n_j);
    epsilon_ij  = rnorm(N); 

# Generate regressors
    k   = dim(corr[[1]])[1]; 
    X_j = rmvnorm(q,rep(0,k),sigma=diag(k));
    X_ij= matrix(0,N,k);
    Z   = matrix(0,N,k);
    ub  =0;

    if (design==3 || design==13){
    for (j in 1:q){
        lb = ub + 1;
        ub = ub + n_j[j];
        X_ij[lb:ub,] = rmvnorm(n_j[j],sigma=corr[[j]]);
        Z[lb:ub,]    = X_ij[lb:ub,] + X_j[rep(j,n_j[j]),];   
    }
    } else if (design==4 || design==14){
    for (j in 1:q){
        lb = ub + 1;
        ub = ub + n_j[j];
        if (j<= q/2){
         Z[lb:ub,] = rmvnorm(n_j[j],mean=c(2,4),sigma=corr[[1]]);
        } else { Z[lb:ub,] = rmvnorm(n_j[j],mean=c(-4,-2),sigma=corr[[2]]);}    
        } 
    }    

# Error term 
if (design==3){
    U    = (Z[,1]^2)*(delta_j + epsilon_ij);
} else if (design==4) {
    U    = (Z[,1]+Z[,2])^2*(delta_j + epsilon_ij);
} else if (design==13){
    U    = delta_j + (epsilon_ij*X_ij[,1]^2+X_ij[,1]^2-1);
} else if (design==14){
    U    = delta_j + (epsilon_ij*(X_ij[,1]+1)^2+(X_ij[,1]+1)^2-1);
}

# Generate Outcome
    beta = cbind(c(beta1,1,1))
    if (design==4 || design==14){beta = cbind(c(beta1,2)); }
    Y    = 1 + Z%*%beta + U;

# Return the data frame
    df  = data.frame(cluster,unit,Y,Z);
    return(df);
}

#-------------------------------------------------------------------
gen.corr <-function(k,q,homo=0){
#-------------------------------------------------------------------
# X-plain: Generates data for the simulations
#-------------------------------------------------------------------
# INPUTS: - homo: 1 (homogeneous) and 0 (heterogenous)
#         - k: dimension of the covariates
#------------------------------------------------------------------
# RETURNS: either a kxk correlation mat or q kxk correlations mat 
#------------------------------------------------------------------
    # homogeneous case
    if (homo==1){
        corr<-rcorr(k);
    }

    # heterogeneous case
    if (homo==0){
        corr<-vector("list", q);
        for (j in 1:q){
           corr[[j]]<-rcorr(k);   
        }
    }

    return(corr)
}

#-------------------------------------------------------------------
gen.Z <-function(q,n,corr,homo=0){
#-------------------------------------------------------------------
# X-plain: Generates data for the simulations
#-------------------------------------------------------------------
# INPUTS: - homo: 1 (homogeneous) and 0 (heterogenous)
#         - corr: correlation matrix from gen.corr() - Determines dim of Z
#         - q: the number of clusters
#         - n: the reference sample size within cluster
#------------------------------------------------------------------
# RETURNS: covariates Z of dimension n times q. 
#------------------------------------------------------------------
    # homogeneous case
    if (homo==1){
        Z<- rmvnorm(n*q,sigma=corr);
    }

    # heterogeneous case
    if (homo==0){
        k   = dim(corr[[1]])[1]; 
        Z_j = rmvnorm(q,rep(0,k),sigma=diag(k));
        Z   = matrix(0,n*q,k);

        for (j in 1:q){
            lb = (j-1)*n+1
            ub = j*n
            Z[lb:ub,] = rmvnorm(n,sigma=corr[[j]]) + Z_j[rep(j,n),];   
        }
    }
    return(Z)
}



#-------------------------------------------------------------------
gen.covariates  <- function(data,fe=1){
#-------------------------------------------------------------------
# X-plain: Computes covariates for the CCE variance estimator 
#-------------------------------------------------------------------
# INPUTS: - data: the entire data frame
#         - fe: 1 (if fixed effects) and 0 (if simple)
#------------------------------------------------------------------
# RETURNS: - A matrix with all covarites
#------------------------------------------------------------------
# Generate useful variables
    clu  =  data$cluster;
    Z    =  data[,c(-(1:3))];
    N    =  length(clu);
    q    =  length(unique(clu));
    ones =  matrix(1,N,1);
    c.dummies = dummy(data$cluster)[,-1];
    if (fe==1){
         covariates= as.matrix(cbind(ones,Z,c.dummies));
    } else {
         covariates= as.matrix(cbind(ones,Z));
    }
return(covariates)
}

#------------------------------------------------------------------
cce.stata <- function(X,reg,cluster,dof="stata") {
# Computes the robust cluster variance estimator implemented in Stata
# ... as on p.16 in NBER WP 18478 by Imbens and Kolesar
#------------------------------------------------------------------
# INPUTS: - X:       N-by-L matrix of regressors
#         - reg:     outcome of a model fit (to obtain residuals and for sandwich)
#         - cluster: N-by-1 vector denoting cluster membership of observations
#                    clusters can be denoted by anything, not just consecutive numbers
#         - dof:    degrees of freedom correction: "stata" default. Also "bch". 
#------------------------------------------------------------------
# WARNING: if constant is included in the regression 'reg', it must be supplied
#          ... as part of X
#------------------------------------------------------------------
# RETURNS: - Estimate of asymptotic var-cov matrix divided by N 
#------------------------------------------------------------------
    N <- length(cluster)             # Number of observations
    S <- length(unique(cluster))     # Number of clusters
    L <- dim(X)[2]                   # Number of regressors
    ## In R, a vector is dimensionless. As such, X[1,] returns an error.
    ## Use the following to introduce consistency with matrix X references.
    if (is.null(L)) {
        L <- 1
        X <- as.matrix(X)
    }
    res <- as.vector(residuals(reg))  # Model residuals

    ## Construct meat matrix
    mymeat <- matrix(rep(0,L^2),L,L)
    for (s in unique(cluster)) {
        tag <- cluster==s
        if (sum(tag)==1) {
            Xs <- t(X[tag,])
        } else {
            Xs <- X[tag,]
        }
        mymeat <- mymeat + crossprod(crossprod(res[tag],Xs))
    }
    mymeat <- mymeat/N
    # DFC as in BCH11 pg. 142. or as on p.16 in NBER WP 18478 by Imbens and Kolesar
    dfc <- (S/(S-1))*(dof=="bch")+(N/(N-L))*(dof=="HC1")+((N-1)/(N-L))*(S/(S-1))*(dof=="stata");

    Vstata <- dfc*sandwich(reg,meat.=mymeat)
    return(Vstata)
}

#------------------------------------------------------------------
Wild.Boot.NS <-function(data,beta.null,alpha,fe=1,Ws="rade"){
#-------------------------------------------------------------------
# X-Plain: Computes the wild bootstrap for Non-Studentized statistic
#          the code only allows for inference on the first covariate.
#-------------------------------------------------------------------
# INPUTS: - data: the original data  
#         - beta.null: the value of beta under H_0
#         - alpha: significance level
#         - fe: 1 (includes clusters fixed effects) and 0 (no fe)
#         - Ws: could be "rade" or "mammen"
#------------------------------------------------------------------
    clu <- data$cluster;
    N   <- length(clu);
    q   <- length(unique(clu));
    X1  <- data$X1;
    
    # get rademacher random variables (could choose mammen too)
    rade <- WB.weights(q,Ws);
    Bs   <- dim(rade)[2];
   
    # formulas for regression with or without fixed effects
    if (fe==1){
        con.form <-Y-beta.null*X1 ~ .-X1-unit-cluster+factor(cluster);
        B.form <- B.Y ~ .-Y-unit-cluster+factor(cluster);
    } else {
        con.form <-Y-beta.null*X1 ~ .-X1-unit-cluster;
        B.form <- B.Y~ .-Y-unit-cluster;
    }

    # Step 1: get constrained estimator and residuals
    con.reg<- lm(con.form,data);
    resi   <- residuals(con.reg);
    Y.hat  <- fitted.values(con.reg)+beta.null*X1;

    B.beta<- matrix(0,Bs);

    # Step 2: do the bootstrap using rademacher weights
    for (s in 1:Bs){
        rademacher<- cbind(rep(c(rade[,s]),c(table(clu))));
        B.U     <- resi*rademacher;
        # Compute the bootstrapped Y's
        B.Y <- Y.hat + B.U;
        # Run LS using the bootstrapped data (Y.ast,original cov).
        B.reg <-lm(B.form,data=data);
        B.beta[s]<-coef(B.reg)[2];
    }

    # Step 3: compute quantiles of bootstrapped test statistic
    B.stat   <- sqrt(N)*(B.beta-beta.null);
    B.abs    <- sqrt(N)*abs(B.beta-beta.null);

    cH.wild  <- quantile(B.stat,1-alpha/2,type=1);
    cL.wild  <- quantile(B.stat,alpha/2,type=1);
    cv.wild  <- quantile(B.abs,1-alpha,type=1);

    return(list(cH=cH.wild,cL=cL.wild,cv=cv.wild))
    #return(cv.wild)
}

#------------------------------------------------------------------
Wild.Boot.S <-function(data,beta.null,alpha,fe=1,Ws="rade"){
#-------------------------------------------------------------------
# X-Plain: Computes the wild bootstrap for Studentized statistic
#          the code only allows for inference on the first covariate.
#-------------------------------------------------------------------
# INPUTS: - data: the original data  
#         - beta.null: the value of beta under H_0
#         - alpha: significance level
#         - fe: 1 (includes clusters fixed effects) and 0 (no fe)
#         - Ws: could be "rade" or "mammen"
#------------------------------------------------------------------
    clu   <- data$cluster;
    N     <- length(clu);
    q     <- length(unique(clu));
    X1    <- data$X1;
    ones  <- matrix(1,N,1);
    all.Z  <- data[,-c(1:3)];
    c.dummies <-dummy(data$cluster)[,-1];
    c.alpha <- alpha;
    if (2^(1-q)<alpha){
    c.alpha  <- alpha - 2^(1-q);
    }

    # get Wild Bootstrap random variables
    weights<-WB.weights(q,Ws);
    Bs     <- dim(weights)[2];
 
    # formulas for regression with or without fixed effects
    if (fe==1){
        con.form <-Y-beta.null*X1 ~ .-X1-unit-cluster+factor(cluster);
        B.form <- B.Y~ .-Y-unit-cluster+factor(cluster);
        covariates<-as.matrix(cbind(ones,all.Z,c.dummies));
    } else {
        con.form <-Y-beta.null*X1 ~ .-X1-unit-cluster;
        B.form <- B.Y~ .-Y-unit-cluster;
        covariates<-as.matrix(cbind(ones,all.Z));
    }

    # Step 1: get constrained estimator and residuals
    con.reg <- lm(con.form,data);
    resi    <- residuals(con.reg);
    Y.hat   <- fitted.values(con.reg)+beta.null*X1;

    B.beta<- matrix(0,Bs); 
    B.cce  <-matrix(0,Bs); 

    # Step 2: do the bootstrap using rademacher weights
    for (s in 1:Bs){
        rademacher<- cbind(rep(c(weights[,s]),c(table(clu))));
        B.U     <- resi*rademacher;
        # Compute the bootstrapped Y's
        B.Y <- Y.hat + B.U;
        # Run LS using the bootstrapped data (Y.ast,original cov).
        B.reg      <-lm(B.form,data=data);
        B.beta[s]  <-coef(B.reg)[2];
        B.cce[s]  <- sqrt(N*cce.stata(covariates,B.reg,clu,"stata")[2,2]);
    }

    # Step 3: compute quantiles of bootstrapped test statistic

    B.stat   <- sqrt(N)*(B.beta-beta.null)/B.cce;
    B.abs    <- sqrt(N)*abs(B.beta-beta.null)/B.cce;

    cH.wild  <- quantile(B.stat,1-alpha/2,type=1);
    cL.wild  <- quantile(B.stat,alpha/2,type=1);
    cv.wild  <- quantile(B.abs,1-alpha,type=1);

    cH.correct  <- quantile(B.stat,1-c.alpha/2,type=1);
    cL.correct  <- quantile(B.stat,c.alpha/2,type=1);
    cv.correct  <- quantile(B.abs,1-c.alpha,type=1);

    return(list(cH=cH.wild,cL=cL.wild,cv=cv.wild,co.H=cH.correct,co.L=cL.correct,co.cv=cv.correct))
}

#------------------------------------------------------------------
Wild.Boot <-function(data,beta.null,alpha,fe=1,Ws="rade",method=0){
#-------------------------------------------------------------------
# X-Plain: Computes the wild bootstrap for Studentized and Non-Studentized statistic
#          This code is for ONE covariate Z in the DGP 
#-------------------------------------------------------------------
# INPUTS: - data: the original data  
#         - beta.null: the value of beta under H_0
#         - alpha: significance level
#         - fe: 1 (includes clusters fixed effects) and 0 (no fe)
#         - Ws: could be "rade" or "mammen"
#------------------------------------------------------------------
    clu     = data$cluster;
    N       = length(clu);
    q       = length(unique(clu));
    Z       = data$Z;
    ones    = matrix(1,N,1);
    all.Z   = Z;
    c.dummies= dummy(data$cluster)[,-1];
    c.alpha = alpha;
    
    if (2^(1-q)<alpha){
        c.alpha = alpha - 2^(1-q);
    }
    
    # get rademacher random variables (could choose mammen too)
    rade = WB.weights(q,Ws,method);
    Bs   = dim(rade)[2];
   
    # formulas for regression with or without fixed effects
    if (fe==1){
        con.form  = Y-beta.null*Z ~ .-unit-cluster+factor(cluster);
        B.form    = B.Y ~ .-Y-unit-cluster+factor(cluster);
        covariates= as.matrix(cbind(ones,all.Z,c.dummies));
    } else {
        con.form  = Y-beta.null*Z ~ .-unit-cluster;
        B.form    = B.Y~ .-Y-unit-cluster;
        covariates=as.matrix(cbind(ones,all.Z));
    }

    # Step 1: get constrained estimator and residuals
    con.reg = lm(con.form,data);
    resi    = residuals(con.reg);
    Y.hat   = fitted.values(con.reg)+beta.null*Z;

    B.beta  = matrix(0,Bs);
    B.cce   = matrix(0,Bs); 

    # Step 2: do the bootstrap using rademacher/Mammen weights
    for (s in 1:Bs){
        rademacher = cbind(rep(c(rade[,s]),c(table(clu))));
        B.U        = resi*rademacher;
        # Compute the bootstrapped Y's
        B.Y  = Y.hat + B.U;
        # Run LS using the bootstrapped data (Y.ast,original cov).
        B.reg       = lm(B.form,data=data);
        B.beta[s]   = coef(B.reg)[2];
        B.cce[s]    = sqrt(N*cce.stata(covariates,B.reg,clu,"stata")[2,2]);
    }

    # Step 3: compute quantiles of bootstrapped test statistic
    # NON-STUDENTIZED
    NS.stat  = sqrt(N)*(B.beta-beta.null);
    NS.abs   = sqrt(N)*abs(B.beta-beta.null);

    cH.NS    = quantile(NS.stat,1-alpha/2,type=1);
    cL.NS    = quantile(NS.stat,alpha/2,type=1);
    cv.NS    = quantile(NS.abs,1-alpha,type=1);

    # STUDENTIZED
    S.stat  = sqrt(N)*(B.beta-beta.null)/B.cce;
    S.abs   = sqrt(N)*abs(B.beta-beta.null)/B.cce;

    cH.S    = quantile(S.stat,1-alpha/2,type=1);
    cL.S    = quantile(S.stat,alpha/2,type=1);
    cv.S    = quantile(S.abs,1-alpha,type=1);

    # STUDENTIZED - CORRECTED
    cH.Sco  = quantile(S.stat,1-c.alpha/2,type=1);
    cL.Sco  = quantile(S.stat,c.alpha/2,type=1);
    cv.Sco  = quantile(S.abs,1-c.alpha,type=1);

    return(list(cH.NS=cH.NS,cL.NS=cL.NS,cv.NS=cv.NS,cH.S=cH.S,cL.S=cL.S,cv.S=cv.S,cH.Sco=cH.Sco,cL.Sco=cL.Sco,cv.Sco=cv.Sco,Bs=Bs,cal=c.alpha))
}

#------------------------------------------------------------------
Wild.Boot.k <-function(data,beta.null,alpha,fe=1,Ws="rade",method=0){
#-------------------------------------------------------------------
# X-Plain: Computes the wild bootstrap for Studentized statistic
#          the code only allows for inference on the first covariate 
#          but allows for many other covariates in the regression
#-------------------------------------------------------------------
# INPUTS: - data: the original data  
#         - beta.null: the value of beta under H_0
#         - alpha: significance level
#         - fe: 1 (includes clusters fixed effects) and 0 (no fe)
#         - Ws: could be "rade" or "mammen"
#------------------------------------------------------------------
    clu   = data$cluster;
    N     = length(clu);
    q     = length(unique(clu));
    X1    = data$X1;
    ones  = matrix(1,N,1);
    all.Z = data[,-c(1:3)];
    c.dummies = dummy(data$cluster)[,-1];
    c.alpha   = alpha;
    if (2^(1-q)<alpha){
    c.alpha   = alpha - 2^(1-q);
    }

    # get Wild Bootstrap random variables
    weights=WB.weights(q,Ws,method);
    Bs     = dim(weights)[2];
 
    # formulas for regression with or without fixed effects
    if (fe==1){
        con.form  = Y-beta.null*X1 ~ .-X1-unit-cluster+factor(cluster);
        B.form    = B.Y~ .-Y-unit-cluster+factor(cluster);
        covariates= as.matrix(cbind(ones,all.Z,c.dummies));
    } else {
        con.form  = Y-beta.null*X1 ~ .-X1-unit-cluster;
        B.form    = B.Y~ .-Y-unit-cluster;
        covariates=as.matrix(cbind(ones,all.Z));
    }

    # Step 1: get constrained estimator and residuals
    con.reg = lm(con.form,data);
    resi    = residuals(con.reg);
    Y.hat   = fitted.values(con.reg)+beta.null*X1;

    B.beta = matrix(0,Bs); 
    B.cce  = matrix(0,Bs); 

    # Step 2: do the bootstrap using rademacher weights
    for (s in 1:Bs){
        rademacher = cbind(rep(c(weights[,s]),c(table(clu))));
        B.U        = resi*rademacher;
        # Compute the bootstrapped Y's
        B.Y = Y.hat + B.U;
        # Run LS using the bootstrapped data (Y.ast,original cov).
        B.reg     = lm(B.form,data=data);
        B.beta[s] = coef(B.reg)[2];
        B.cce[s]  = sqrt(N*cce.stata(covariates,B.reg,clu,"stata")[2,2]);
    }

    # Step 3: compute quantiles of bootstrapped test statistic
    # NON-STUDENTIZED
    NS.stat  = sqrt(N)*(B.beta-beta.null);
    NS.abs   = sqrt(N)*abs(B.beta-beta.null);

    cH.NS    = quantile(NS.stat,1-alpha/2,type=1);
    cL.NS    = quantile(NS.stat,alpha/2,type=1);
    cv.NS    = quantile(NS.abs,1-alpha,type=1);

    # STUDENTIZED
    S.stat  = sqrt(N)*(B.beta-beta.null)/B.cce;
    S.abs   = sqrt(N)*abs(B.beta-beta.null)/B.cce;

    cH.S    = quantile(S.stat,1-alpha/2,type=1);
    cL.S    = quantile(S.stat,alpha/2,type=1);
    cv.S    = quantile(S.abs,1-alpha,type=1);

    # STUDENTIZED - CORRECTED
    cH.Sco  = quantile(S.stat,1-c.alpha/2,type=1);
    cL.Sco  = quantile(S.stat,c.alpha/2,type=1);
    cv.Sco  = quantile(S.abs,1-c.alpha,type=1);

    return(list(cH.NS=cH.NS,cL.NS=cL.NS,cv.NS=cv.NS,cH.S=cH.S,cL.S=cL.S,cv.S=cv.S,cH.Sco=cH.Sco,cL.Sco=cL.Sco,cv.Sco=cv.Sco,Bs=Bs,cal=c.alpha))

}

#-------------------------------------------------------------------
WB.weights  <- function(q,Ws="rade",method=1){
#-------------------------------------------------------------------
# X-plain: Computes weights for the wild bootstrap 
#-------------------------------------------------------------------
# INPUTS: - q: the dimension of the data (q point estimators)
#         - Ws: could be "rade" or "mammen"
#         - method: 1 (stochastic approx) or 0 (exact)
#           note, exact only available for q<11. 
#------------------------------------------------------------------
# RETURNS: - G: a B \times q matrix.
#------------------------------------------------------------------
Bs=500*(q<10) + 1000*(q>=10);

if (Ws=="mammen"){
    u <-runif(q*Bs);
    c <-(sqrt(5)+1)/(2*sqrt(5));
    weights <- matrix((u<=c)*(1-sqrt(5))/2+(u>c)*(1+sqrt(5))/2,q,Bs);
}

if (Ws=="rade"){
   if (q<11){
        if (method==0){
        e   <- c(1,-1); Bs<-2^q;
        weights<-t(expand.grid(rep(list(e),q)));
        } else {
            weights <- matrix(2*((runif(q*Bs)>1/2)-1/2),q,Bs);
            }
    } else {
        weights <- matrix(2*((runif(q*Bs)>1/2)-1/2),q,Bs);     
    }
}

return(weights)
}
#-------------------------------------------------------------------
gaussmix <- function(nsim,w=0.9,mean1=-1/9,mean2=1,sd1=1,sd2=2){
#-------------------------------------------------------------------
# X-plain Simulates data from a mixture of two normals 
#-------------------------------------------------------------------
# INPUTS: - nsim: number of data to be simulated
#------------------------------------------------------------------
# RETURNS: - y: the simulated data
#------------------------------------------------------------------
   U <- runif(nsim)
   I <- as.numeric(U<w);
   y <- I*rnorm(nsim,mean=mean1,sd=sd1)+
       (1-I)*rnorm(nsim,mean=mean2,sd=sd2);
   return(y)
}
#------------------------------------------------------------------
# END OF FILE
#------------------------------------------------------------------
